#include <iostream>
#include<sys/time.h>
#include<pthread.h>
#include<semaphore.h>
#include<arm_neon.h>
#define N 100
#define C 100
#define P 7
//pth编号从零开始依次为：动态线程、静态线程行分配、静态线程列分配
using namespace std;
float A[N][N];
typedef struct {
    int k;//消去的轮次
    int t_id;//线程id
}threadParam_t;
sem_t sem_main;
sem_t sem_workerstart[P];
sem_t sem_workerend[P];
void* threadFunc(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int k = p->k;//消去的轮次
    int t_id = p->t_id;//线程编号
    int i = k + t_id +1;//获取自己的计算任务
    for(int j = k + 1; j < N; j++) {
        A[i][j] = A[i][j] - A[i][k] * A[k][j];
    }
    A[i][k] = 0;
    pthread_exit(nullptr);
}
void* threadFunc1(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;
    for(int k = 0; k < N; k++) {
        sem_wait(&sem_workerstart[t_id]);//阻塞，等待主线程完成除法操作
        //循环划分任务
        for(int i = k + t_id + 1; i < N; i += P) {
            //消去
            for(int j = k + 1; j < N; j++) {
                A[i][j] = A[i][j] - A[i][k]*A[k][j];
            }
            A[i][k] = 0.0;
        }
        sem_post(&sem_main);//唤醒主线程
        sem_wait(&sem_workerend[t_id]);//阻塞，等待主线程唤醒进入下一轮
    }
    pthread_exit(nullptr);
}
void* threadFunc2(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;
    for(int k = 0; k < N; k++) {
        sem_wait(&sem_workerstart[t_id]);//阻塞，等待主线程完成除法操作
        //循环划分任务
        for(int i = k + 1; i < N; i++) {
            //消去
            for(int j = k + t_id + 1; j < N; j += P) {
                A[i][j] = A[i][j] - A[i][k]*A[k][j];
            }
            A[i][k] = 0.0;
        }
        sem_post(&sem_main);//唤醒主线程
        sem_wait(&sem_workerend[t_id]);//阻塞，等待主线程唤醒进入下一轮
    }
     pthread_exit(nullptr);
}
void display() {
    for(int i = 0; i < N; i++) {
        for(int j = 0; j < N; j++) {
            cout<<A[i][j]<<" ";
        }
        cout<<endl;
    }
}//check outcome
void init() {
    //cout<<"Initializing...";
    //cout<<"Initializing ...\n";
    for(int i = 0; i < N; i++) {
        A[i][i] = 1.0;
    }
    for(int r = 0; r < 5 * N; r++) {
        int i1 = rand()%N;
        int i2 = rand()%N;
        float rate = rand()%10/10.0;;
        if(i1 != i2) {
            for(int j = 0; j < N; j++) {
                A[i1][j] += rate*A[i2][j];
            }
        }
    }
    //cout<<"Finished. \n";
}
//平凡算法
void normal() {
    struct timeval h, t;
    float time = 0.0;
    //Normal
    cout<<"Normal:(N = "<<N<<")\n";
    for(int r = 0; r < C; r++) {
        init();
        gettimeofday(&h, NULL);
        for(int k = 0; k < N; k++) {
            for(int j = k + 1; j < N; j++) {
                A[k][j] = A[k][j]/A[k][k];
            }
            A[k][k] = 1.0;
            for(int i = k + 1; i < N; i++) {
                for(int j = k + 1; j < N; j++) {
                    A[i][j] = A[i][j] - A[i][k]*A[k][j];
                }
                A[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME: "<<time/C<<" ms"<<endl;
}
//动态线程
void pth() {
    struct timeval h, t;
    float time = 0.0;
    cout<<"Pthread:(N = "<<N<<")\n";
    for(int r = 0; r < C; r++) {
        init();
        gettimeofday(&h, NULL);
        for(int k = 0; k < N; k++) {
            //主线程做除法操作
            for(int j = k + 1; j < N; j++) {
                A[k][j] = A[k][j]/A[k][k];
            }
            A[k][k] = 1.0;
            //创建工作线程，进行消去操作
            int worker_count = N - 1 - k;
            pthread_t* handles = (pthread_t*)malloc(N*sizeof(pthread_t));
            threadParam_t* param = (threadParam_t*)malloc(N*sizeof(threadParam_t));
            //分配任务
            for(int t_id = 0; t_id < worker_count; t_id++) {
                param[t_id].k = k;
                param[t_id].t_id = t_id;
            }
            //创建线程
            for(int t_id = 0; t_id <worker_count; t_id++) {
                pthread_create(handles + t_id, nullptr, threadFunc, param + t_id);
            }
            //主线程挂起等待所有工作线程完成此轮消去工作
            for(int t_id = 0; t_id < worker_count; t_id++) {
                pthread_join(handles[t_id], nullptr);
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME: "<<time/C<<" ms"<<endl;
}
//静态线程行分配
void pth1() {
    struct timeval h, t;
    float time = 0.0;
    cout<<"Pthread1:(N = "<<N<<" and P = "<<P<<")\n";
    for(int r = 0; r < C; r++) {
        init();
        gettimeofday(&h, NULL);
        //初始化信号量
        sem_init(&sem_main, 0, 0);
        for(int i = 0; i < P; i++) {
            sem_init(&sem_workerstart[i], 0, 0);
            sem_init(&sem_workerend[i], 0, 0);
        }
        //创建线程
        pthread_t handles[P];
        threadParam_t param[P];
        for(int t_id = 0; t_id < P; t_id++) {
            param[t_id].t_id = t_id;
            pthread_create(handles + t_id, nullptr, threadFunc1, param + t_id);
        }
        for(int k = 0; k < N; k++) {
            //主线程做除法操作
            for(int j = k + 1; j < N; j++) {
                A[k][j] = A[k][j]/A[k][k];
            }
            A[k][k] = 1.0;
            //开始唤醒工作线程
            for(int t_id = 0; t_id < P; t_id++) {
                sem_post(&sem_workerstart[t_id]);
            }
            //主线程睡眠
            for(int t_id = 0; t_id < P; t_id++) {
                sem_wait(&sem_main);
            }
            //主线程再次唤醒工作线程
            for(int t_id = 0; t_id < P; t_id++) {
                sem_post(&sem_workerend[t_id]);
            }
        }
        for(int t_id = 0; t_id < P; t_id++) {
            pthread_join(handles[t_id], nullptr);
        }
        sem_destroy(&sem_main);
        for(int i = 0; i < P; i++) {
            sem_destroy(&sem_workerstart[i]);
            sem_destroy(&sem_workerend[i]);
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME: "<<time/C<<" ms"<<endl;
}
//静态线程列分配
void pth2() {
    struct timeval h, t;
    float time = 0.0;
    cout<<"Pthread2:(N = "<<N<<" and P = "<<P<<")\n";
    for(int r = 0; r < C; r++) {
        init();
        gettimeofday(&h, NULL);
        //初始化信号量
        sem_init(&sem_main, 0, 0);
        for(int i = 0; i < P; i++) {
            sem_init(&sem_workerstart[i], 0, 0);
            sem_init(&sem_workerend[i], 0, 0);
        }
        //创建线程
        pthread_t handles[P];
        threadParam_t param[P];
        for(int t_id = 0; t_id < P; t_id++) {
            param[t_id].t_id = t_id;
            pthread_create(handles + t_id, nullptr, threadFunc2, param + t_id);
        }
        for(int k = 0; k < N; k++) {
            //主线程做除法操作
            for(int j = k + 1; j < N; j++) {
                A[k][j] = A[k][j]/A[k][k];
            }
            A[k][k] = 1.0;
            //开始唤醒工作线程
            for(int t_id = 0; t_id < P; t_id++) {
                sem_post(&sem_workerstart[t_id]);
            }
            //主线程睡眠
            for(int t_id = 0; t_id < P; t_id++) {
                sem_wait(&sem_main);
            }
            //主线程再次唤醒工作线程
            for(int t_id = 0; t_id < P; t_id++) {
                sem_post(&sem_workerend[t_id]);
            }
        }
        for(int t_id = 0; t_id < P; t_id++) {
            pthread_join(handles[t_id], nullptr);
        }
        sem_destroy(&sem_main);
        for(int i = 0; i < P; i++) {
            sem_destroy(&sem_workerstart[i]);
            sem_destroy(&sem_workerend[i]);
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME: "<<time/C<<" ms"<<endl;
}
int main() {
    normal();
    //display();
    pth();
    //display();
    pth1();
    //display();
    pth2();
    //display();
}